#####################
## Distance to road function
#####################


dist2road <- function(mod.vil){
X <- matrix(NA, ncol=ncol(mod.vil$model[,-1]),nrow=1000)
medians <- apply(mod.vil$model[,-1], MARGIN=2,FUN=median)
counter <- 1:ncol(X)
for(i in counter){
X[,i] <- as.numeric(medians)[i]
}
X[,1] <- seq(0,max(mod.vil$model[,2]),length=1000)
X <- as.data.frame(X)
#X[,which(colnames(mod.vil$model[,-1])%in%"OBLAST_NAM")] <- as.data.frame(medians)["OBLAST_NAM",]
#X[,which(colnames(mod.vil$model[,-1])%in%"OBLAST_NAM")] <- unique(data$OBLAST_NAM)[7]
colnames(X) <- colnames(mod.vil$model[,-1])
pred <- predict.gam(mod.vil,newdata=X,type="link",se.fit=T)
Y <- 1/(1+exp(-pred$fit))
Y.u <- 1/(1+exp(-(pred$fit+1.96*pred$se.fit)))
Y.l <- 1/(1+exp(-(pred$fit-1.96*pred$se.fit)))
return(list(X=X,Y=Y,Y.l=Y.l,Y.u=Y.u))
}

#####################
## Distance to road function (conditional)
#####################



dist2road.cond <- function(mod.vil){
X <- matrix(NA, ncol=ncol(mod.vil$model[,-1]),nrow=1000)
medians <- apply(mod.vil$model[,-1], MARGIN=2,FUN=median)
inter <- grep("R_",colnames(mod.vil$model[,-1]))
pre.inter <- gsub("R_","",colnames(mod.vil$model[,-1])[inter])

counter <- 1:ncol(X)
for(i in counter){
X[,i] <- as.numeric(medians)[i]
}
X[,1] <- seq(0,max(mod.vil$model[,2]),length=1000)
X <- as.data.frame(X)
colnames(X) <- colnames(mod.vil$model[,-1])

# No attack
pred <- predict.gam(mod.vil,newdata=X,type="link",se.fit=T)
Y0 <- 1/(1+exp(-pred$fit))
Y0.u <- 1/(1+exp(-(pred$fit+1.96*pred$se.fit)))
Y0.l <- 1/(1+exp(-(pred$fit-1.96*pred$se.fit)))

# Attack
X$REBEL.b.t1 <- 1
X[,inter] <- X[,pre.inter]

pred <- predict.gam(mod.vil,newdata=X,type="link",se.fit=T)
Y1 <- 1/(1+exp(-pred$fit))
Y1.u <- 1/(1+exp(-(pred$fit+1.96*pred$se.fit)))
Y1.l <- 1/(1+exp(-(pred$fit-1.96*pred$se.fit)))
return(list(X=X,Y0=Y0,Y0.l=Y0.l,Y0.u=Y0.u,Y1=Y1,Y1.l=Y1.l,Y1.u=Y1.u))
}



######################
## Closest event 
######################

closest <- function(y,W){
	dist <- c()
	W <- as.matrix(W,nrow=nrow(W),ncol=ncol(W))
	diag(W) <- max(W)
	dist <- ifelse(rep(sum(y)>0,nrow(W)),apply(cbind(W[,y==1],max(W)),1,min),rep(max(W),nrow(W)))
	return(dist)}
	
	
######################
## Labels
######################

.simpleCap <- function(x){
    s <- strsplit(x, " ")[[1]]
    paste(toupper(substring(s, 1,1)), substring(s, 2), sep="", collapse=" ")}


######################
## Round breaks for color palette
######################

round.breaks <- function(cols){
	breaks <- names(attr(cols, "table"))
	breaks <- strsplit(x=breaks, c(","))
	breaks <- unlist(breaks)
	breaks <- sub(pattern="[", replacement="", x=breaks,fixed=T)
	breaks <- sub(pattern="]", replacement="", x=breaks,fixed=T)
	breaks <- sub(pattern="(", replacement="", x=breaks,fixed=T)
	breaks <- sub(pattern=")", replacement="", x=breaks,fixed=T)
	n <- length(breaks)
	breaks <- round(as.numeric(breaks),4)
	odd <- subset(subset(0:(n-2),0:(n-2)%%2==1),is.na(subset(0:(n-2),0:(n-2)%%2==1))==0)
	even <- subset(subset(2:(n-2),0:(n-2)%%2==0),is.na(subset(2:(n-2),0:(n-2)%%2==0))==0)
	breaks <- c(paste("[",breaks[odd]," , ",breaks[even],")",sep=""),paste("[",breaks[n-1]," , ",breaks[n],"]",sep=""))
	return(breaks) }

######################
## Degree to kilometer
######################

d2km <- function(d){
out <- d*60*1.852
return(out)
}


######################
## Kilometer to degree
######################

km2d <- function(km){
out <- (km/1.852)/60
return(out)
}


######################
## Neat scalebar
######################

scalebar <- function(loc,length,unit="km",division.cex=.7,...) {
if(missing(loc)) stop("loc is missing")
if(missing(length)) stop("length is missing")
x <- c(0,length/c(4,2,4/3,1),length*1.1)+loc[1]
y <- c(0,length/(10*3:1))+loc[2]
cols <- rep(c("black","white"),2)
for (i in 1:4) rect(x[i],y[1],x[i+1],y[2],col=cols[i])
for (i in 1:5) segments(x[i],y[2],x[i],y[3])
labels <- x[c(1,3)]-loc[1]
labels <- append(round(d2km(labels),0),paste(round(d2km(x[5]-loc[1]),0),unit))
text(x[c(1,3,5)],y[4],labels=labels,adj=.5,cex=division.cex,offset=.1,pos=3)
}

######################
## AUC in-sample (GAM)
######################

AUC <- function(y, fitted){
    	n1 <- sum(y)
    	n <- length(y)
    	out <- (mean(rank(fitted)[y == 1]) - (n1 + 1)/2)/(n - n1)
    	return(out)}


######################
## AUC out-sample (GAM)
######################

AUC.out <- function(mod,data.out){
	mat.X <- data.out[,names(mod$model)]
	pred.out <- predict.gam(mod,newdata=mat.X[,-1],type="response")
	out <- AUC(y=data.out[,names(mod$model)[1]], fitted=c(pred.out))
	return(out)}


######################
## Simulation map
######################

simMap <- function(mod,data,W,city,title,box.col="blue",window=.5){
	
## Extract cross-section
X <- subset(data,data$YRMO==200812)

## Village names
NAMES <- c()
for(i in 1:length(X$NAME)){
NAMES[i] <- .simpleCap(as.character(X$NAME[i]))}

## Baseline: No violence
set.seed(12345)
X$REBEL.b.t1 <- 0
X <- X[,names(mod$model[,-1])]
if(sum(grepl("D.",names(X),fixed=T))>1){X[,1] <- closest(X$REBEL.b.t1,W=W)}
pred0 <- predict.gam(mod,newdata=X,type="response",se.fit=T)

## Counterfactual: Violence in village
X <- subset(data,data$YRMO==200812)
set.seed(12345)
X$REBEL.b.t1 <- 0
v1 <- city
X$REBEL.b.t1[c(v1)] <- 1
X <- X[,names(mod$model[,-1])]
if(sum(grepl("D.",names(X),fixed=T))>1){X[,1] <- closest(X$REBEL.b.t1,W=W)}
pred1 <- predict.gam(mod,newdata=X,type="response",se.fit=T)

## Difference
ATT <- pred1$fit-pred0$fit
var <- as.numeric(ATT)
classes_fx <- classIntervals(c(var,-6:-1), n=6, style = "fixed", fixedBreaks=c(0,.01, .05, .1,.5, 1),rtimes = 1)


## Colors
#br.palette <- colorRampPalette(c("white", "darkorange","red","brown"), space = "rgb")
br.palette <- colorRampPalette(c("white", "blue","darkblue"), space = "rgb")
pal <- br.palette(5)
cols <- findColours(classes_fx, pal)
breaks <- round.breaks(cols)

## Big plot
X. <- X
X. <- X[order(var),]
cols. <- cols[order(var)]
var. <- var[order(var)]

par(mar=c(.5,.05,.05,10),xpd=T)
plot(x=X.$LONG[var.>.01],y=X.$LAT[var.>.01],col=cols.[var.>.01],cex=1.1,pch=16,axes=F,xlim=c(min(X$LONG),max(X$LONG)),ylim=c(min(X$LAT),max(X$LAT)),xlab="",ylab="")
plot(road,add=T,col="chocolate1",lwd=.3)
points(x=X$LONG[X$REBEL.b.t1==1],y=X$LAT[X$REBEL.b.t1==1],col="gold",cex=1.3)
points(x=X$LONG[X$REBEL.b.t1==1],y=X$LAT[X$REBEL.b.t1==1],col="black",cex=1.5)
plot(map,add=T,border="black",lwd=.5)
par(xpd=NA)
box(lty = "solid", lwd=1, col = box.col,which="figure")
legend(x=48,y=45,cex=.9,fill=attr(cols,"palette"),bty="n",legend=paste(breaks,"   N = ",as.numeric(attr(cols, "table")),sep=""),title="Change in risk of violence",ncol=1)
#text(x=min(X$LONG)-.5,y=min(X$LAT)+.7,labels=bquote(bold(.(title))),cex=1.3,pos=4)
#legend(x="bottomleft",cex=1.1,pch=1,col="red",legend="Outbreak location (t - 1)",bty="n")
legend(x="bottomleft",legend=bquote(bold(.(title))),cex=1.3,bty="n")

segments(x0=X$LONG[v1]-window,x1=X$LONG[v1]-window,y0=X$LAT[v1]+window,y1=X$LAT[v1]-window,col=box.col,lwd=1.5)
segments(x0=X$LONG[v1]+window,x1=X$LONG[v1]+window,y0=X$LAT[v1]+window,y1=X$LAT[v1]-window,col=box.col,lwd=1.5)
segments(x0=X$LONG[v1]-window,x1=X$LONG[v1]+window,y0=X$LAT[v1]-window,y1=X$LAT[v1]-window,col=box.col,lwd=1.5)
segments(x0=X$LONG[v1]-window,x1=X$LONG[v1]+window,y0=X$LAT[v1]+window,y1=X$LAT[v1]+window,col=box.col,lwd=1.5)

# Little plot
xlim <- X$LONG[v1]+c(-window-.25,window-.25)
ylim <- X$LAT[v1]+c(-window,window)
X. <- X[X$LONG>=X$LONG[v1]-window & X$LONG <= X$LONG[v1]+window & X$LAT>=X$LAT[v1]-window & X$LAT <= X$LAT[v1]+window,]
var. <- var[X$LONG>=X$LONG[v1]-window & X$LONG <= X$LONG[v1]+window & X$LAT>=X$LAT[v1]-window & X$LAT <= X$LAT[v1]+window]
cols. <- cols[X$LONG>=X$LONG[v1]-window & X$LONG <= X$LONG[v1]+window & X$LAT>=X$LAT[v1]-window & X$LAT <= X$LAT[v1]+window]
X. <- X.[order(var.),]
cols. <- cols.[order(var.)]
var. <- var.[order(var.)]

par(mar=c(.1,.1,.1,.1),xpd=F)
image(elevation,col=terrain.colors(10),axes=F,xlim=xlim,ylim=ylim)
#plot(x=X.$LONG,y=X.$LAT,col="lightgrey",cex=1,pch=1,axes=F,xlim=c(X$LONG[v1]-window,X$LONG[v1]+window),ylim=c(X$LAT[v1]-window,X$LAT[v1]+window))
points(x=X.$LONG,y=X.$LAT,col="lightgrey",cex=1,pch=1)
plot(road,add=T,col="chocolate3",lwd=1.2)
points(x=X.$LONG[var.>.01],y=X.$LAT[var.>.01],cex=1.2,col=cols.[var.>.01],pch=16)
points(x=X$LONG[X$REBEL.b.t1==1],y=X$LAT[X$REBEL.b.t1==1],col="gold",cex=1.5)
points(x=X$LONG[X$REBEL.b.t1==1],y=X$LAT[X$REBEL.b.t1==1],col="black",cex=1.7)
plot(map,add=T,border="black",lwd=.5)
box(lty = "solid", lwd=2, col = box.col)
#text(x=X$LONG[v1],y=X$LAT[v1],labels=NAMES[v1],cex=1,pos=1)
text(x=X$LONG[v1],y=X$LAT[v1],labels="Grozny",cex=1,pos=1)
}



######################
## Simulation map - conditional
######################


simMap.cond <- function(mod,data,W,city,title,box.col="blue",window=.5){
	
## Indices of interaction terms
inter <- grep("R_",colnames(mod$model[,-1]))
pre.inter <- gsub("R_","",colnames(mod$model[,-1])[inter])
	
## Extract cross-section
X <- subset(data,data$YRMO==200812)

## Village names
NAMES <- c()
for(i in 1:length(X$NAME)){
NAMES[i] <- .simpleCap(as.character(X$NAME[i]))}

## Baseline: No violence
set.seed(12345)
X$REBEL.b.t1 <- 0
X <- X[,names(mod$model[,-1])]
if(sum(grepl("D.",names(X),fixed=T))>1){X[,1] <- closest(X$REBEL.b.t1,W=W)}
X[,inter] <- X[,pre.inter]*X$REBEL.b.t1
pred0 <- predict.gam(mod,newdata=X,type="response",se.fit=T)

## Counterfactual: Violence in village
X <- subset(data,data$YRMO==200812)
set.seed(12345)
X$REBEL.b.t1 <- 0
v1 <- which(NAMES%in%city)
X$REBEL.b.t1[c(v1)] <- 1
X <- X[,names(mod$model[,-1])]
if(sum(grepl("D.",names(X),fixed=T))>1){X[,1] <- closest(X$REBEL.b.t1,W=W)}
X[,inter] <- X[,pre.inter]*X$REBEL.b.t1
pred1 <- predict.gam(mod,newdata=X,type="response",se.fit=T)

## Difference
ATT <- pred1$fit-pred0$fit
var <- as.numeric(ATT)
classes_fx <- classIntervals(c(var,-6:-1), n=6, style = "fixed", fixedBreaks=c(0,.01, .05, .1,.5, 1),rtimes = 1)


## Colors
#br.palette <- colorRampPalette(c("white", "darkorange","red","brown"), space = "rgb")
br.palette <- colorRampPalette(c("white", "blue","darkblue"), space = "rgb")
pal <- br.palette(5)
cols <- findColours(classes_fx, pal)
breaks <- round.breaks(cols)

## Big plot
X. <- X
X. <- X[order(var),]
cols. <- cols[order(var)]
var. <- var[order(var)]

par(mar=c(.5,.05,.05,10),xpd=T)
plot(x=X.$LONG[var.>.01],y=X.$LAT[var.>.01],col=cols.[var.>.01],cex=1.1,pch=16,axes=F,xlim=c(min(X$LONG),max(X$LONG)),ylim=c(min(X$LAT),max(X$LAT)),xlab="",ylab="")
plot(road,add=T,col="chocolate1",lwd=.3)
points(x=X$LONG[X$REBEL.b.t1==1],y=X$LAT[X$REBEL.b.t1==1],col="gold",cex=1.3)
points(x=X$LONG[X$REBEL.b.t1==1],y=X$LAT[X$REBEL.b.t1==1],col="black",cex=1.5)
plot(map,add=T,border="black",lwd=.5)
par(xpd=NA)
box(lty = "solid", lwd=1, col = box.col,which="figure")
legend(x=48,y=45.2,cex=.9,fill=attr(cols,"palette"),bty="n",legend=paste(breaks,"   N = ",as.numeric(attr(cols, "table")),sep=""),title="Change in risk of violence",ncol=1)
secs <- secEpid(mod=mod,W=W, city=city)
legend(x=49.2,y=42.2,cex=.9,col=NA,fill=NA,border=NA,pch=NA,bty="n",legend=paste(round(secs$mean,2)," (",round(secs$sd,2),")",sep=""),title="Secondary cases (s.d.)",ncol=1)
#text(x=min(X$LONG)-.5,y=min(X$LAT)+.7,labels=bquote(bold(.(title))),cex=1.3,pos=4)
legend(x="bottomleft",legend=bquote(bold(.(title))),cex=1.3,bty="n")
max(X$LONG)
#legend(x="bottomleft",cex=1.1,pch=1,col="red",legend="Outbreak location (t - 1)",bty="n")
segments(x0=X$LONG[v1]-window,x1=X$LONG[v1]-window,y0=X$LAT[v1]+window,y1=X$LAT[v1]-window,col=box.col,lwd=1.5)
segments(x0=X$LONG[v1]+window,x1=X$LONG[v1]+window,y0=X$LAT[v1]+window,y1=X$LAT[v1]-window,col=box.col,lwd=1.5)
segments(x0=X$LONG[v1]-window,x1=X$LONG[v1]+window,y0=X$LAT[v1]-window,y1=X$LAT[v1]-window,col=box.col,lwd=1.5)
segments(x0=X$LONG[v1]-window,x1=X$LONG[v1]+window,y0=X$LAT[v1]+window,y1=X$LAT[v1]+window,col=box.col,lwd=1.5)

# Little plot
xlim <- X$LONG[v1]+c(-window,window)
ylim <- X$LAT[v1]+c(-window,window)
X. <- X[X$LONG>=X$LONG[v1]-window & X$LONG <= X$LONG[v1]+window & X$LAT>=X$LAT[v1]-window & X$LAT <= X$LAT[v1]+window,]
var. <- var[X$LONG>=X$LONG[v1]-window & X$LONG <= X$LONG[v1]+window & X$LAT>=X$LAT[v1]-window & X$LAT <= X$LAT[v1]+window]
cols. <- cols[X$LONG>=X$LONG[v1]-window & X$LONG <= X$LONG[v1]+window & X$LAT>=X$LAT[v1]-window & X$LAT <= X$LAT[v1]+window]
X. <- X.[order(var.),]
cols. <- cols.[order(var.)]
var. <- var.[order(var.)]

par(mar=c(.1,.1,.1,.1),xpd=F)
image(elevation,col=terrain.colors(10),axes=F,xlim=xlim,ylim=ylim)
#plot(x=X.$LONG,y=X.$LAT,col="lightgrey",cex=1,pch=1,axes=F,xlim=c(X$LONG[v1]-window,X$LONG[v1]+window),ylim=c(X$LAT[v1]-window,X$LAT[v1]+window))
points(x=X.$LONG,y=X.$LAT,col="lightgrey",cex=1,pch=1)
plot(road,add=T,col="chocolate3",lwd=1.2)
points(x=X.$LONG[var.>.01],y=X.$LAT[var.>.01],cex=1.5,col=cols.[var.>.01],pch=16)
points(x=X$LONG[X$REBEL.b.t1==1],y=X$LAT[X$REBEL.b.t1==1],col="gold",cex=1.5)
points(x=X$LONG[X$REBEL.b.t1==1],y=X$LAT[X$REBEL.b.t1==1],col="black",cex=1.7)
plot(map,add=T,border="black",lwd=.5)
box(lty = "solid", lwd=2, col = box.col)
#text(x=X$LONG[v1],y=X$LAT[v1],labels=NAMES[v1],cex=1,pos=1)
text(x=X$LONG[v1],y=X$LAT[v1],labels="Grozny",cex=1,pos=1)
}






######################
## Simulation map - forecasting
######################


simMap.fcast <- function(mod,data,W,month,title,box.col="blue",window=.5){
	
## Indices of interaction terms
inter <- grep("I_",colnames(mod$model[,-1]))
pre.inter <- gsub("I_","",colnames(mod$model[,-1])[inter])
	
## Extract cross-section
X <- subset(data,data$YRMO==month)

## Village names
NAMES <- c()
for(i in 1:length(X$NAME)){
NAMES[i] <- .simpleCap(as.character(X$NAME[i]))}

## Baseline: No violence
set.seed(12345)
X$REBEL.b.t1 <- 0
X <- X[,names(mod$model[,-1])]
if(sum(grepl("D.",names(X),fixed=T))>1){X[,1] <- closest(X$REBEL.b.t1,W=W)}
X[,inter] <- X[,pre.inter]*X$REBEL.b.t1
pred0 <- predict.gam(mod,newdata=X,type="response",se.fit=T)

## Counterfactual: Violence in village
X <- subset(data,data$YRMO==month)
city <- which(X$REBEL.b.t1==1)
actual <- which(X$REBEL.b==1)
set.seed(12345)
X$REBEL.b.t1 <- 0
v1 <- city
X$REBEL.b.t1[c(v1)] <- 1
X <- X[,names(mod$model[,-1])]
if(sum(grepl("D.",names(X),fixed=T))>1){X[,1] <- closest(X$REBEL.b.t1,W=W)}
X[,inter] <- X[,pre.inter]*X$REBEL.b.t1
pred1 <- predict.gam(mod,newdata=X,type="response",se.fit=T)

## Difference
ATT <- pred1$fit-pred0$fit
var <- as.numeric(ATT)
classes_fx <- classIntervals(c(var,-6:-1), n=6, style = "fixed", fixedBreaks=c(0,.01, .05, .1,.5, 1),rtimes = 1)


## Colors
#br.palette <- colorRampPalette(c("white", "darkorange","red","brown"), space = "rgb")
br.palette <- colorRampPalette(c("white", "blue","darkblue"), space = "rgb")
pal <- br.palette(5)
cols <- findColours(classes_fx, pal)
breaks <- round.breaks(cols)

## Big plot
X. <- X
X. <- X[order(var),]
cols. <- cols[order(var)]
var. <- var[order(var)]

par(mar=c(.5,.05,.05,10),xpd=T)
plot(x=X.$LONG[var.>.01],y=X.$LAT[var.>.01],col=cols.[var.>.01],cex=1.1,pch=16,axes=F,xlim=c(min(X$LONG),max(X$LONG)),ylim=c(min(X$LAT),max(X$LAT)),xlab="",ylab="")
plot(road,add=T,col="chocolate3",lwd=.5)
points(x=X$LONG[X$REBEL.b.t1==1],y=X$LAT[X$REBEL.b.t1==1],col="red",cex=1.5)
points(x=X$LONG[actual],y=X$LAT[actual],col="forestgreen",cex=2)
plot(map,add=T,border="black",lwd=.5)
par(xpd=NA)
box(lty = "solid", lwd=1, col = box.col,which="figure")
legend(x=48,y=45,cex=.9,fill=attr(cols,"palette"),bty="n",legend=paste(breaks,"   N = ",as.numeric(attr(cols, "table")),sep=""),title="Change in risk of violence",ncol=1)
text(x=min(X$LONG)-.5,y=min(X$LAT)+.9,labels=bquote(bold(.(title))),cex=1.3,pos=4)
legend(x="bottomleft",cex=1.1,pch=1,col=c("red","forestgreen"),legend=c("Violence location (t - 1)","Violence location (t)"),bty="n")
legend(x=50,y=min(X$LAT)+.9,legend=c(paste("P(V|V) =",round(mean(var[actual],na.rm=T),4)),paste("P(V|No V) =",round(mean(var[-c(actual)],na.rm=T),4))),bty="n",cex=.8)
}


########################
## Simulation: secondary cases
########################

secEpid <- function(mod,W,city,month=200812,nsims=1000){

## Indices of interaction terms
inter <- grep("R_",colnames(mod$model[,-1]))
pre.inter <- gsub("R_","",colnames(mod$model[,-1])[inter])
	
## Extract cross-section
X <- subset(data,data$YRMO==month)

## Village names
NAMES <- c()
for(i in 1:length(X$NAME)){
NAMES[i] <- .simpleCap(as.character(X$NAME[i]))}
city <- which(NAMES%in%city)

## Baseline: No violence
set.seed(12345)
X$REBEL.b.t1 <- 0
X <- X[,names(mod$model[,-1])]
if(sum(grepl("D.",names(X),fixed=T))>1){X[,1] <- closest(X$REBEL.b.t1,W=W)}
X[,inter] <- X[,pre.inter]*X$REBEL.b.t1
pred0 <- predict.gam(mod,newdata=X,type="response",se.fit=T)

## Counterfactual: Violence in village
X <- subset(data,data$YRMO==month)
X$REBEL.b.t1 <- 0
X$REBEL.b.t1[c(city)] <- 1
X <- X[,names(mod$model[,-1])]
if(sum(grepl("D.",names(X),fixed=T))>1){X[,1] <- closest(X$REBEL.b.t1,W=W)}
X[,inter] <- X[,pre.inter]*X$REBEL.b.t1
pred1 <- predict.gam(mod,newdata=X,type="response",se.fit=T)
attacks <- c()
for(i in 1:nsims){
preds <- rbinom(n=4033,size=1,prob=pred1$fit-pred0$fit)
attacks[i] <- sum(preds[-city],na.rm=T)}
return(list(mean=mean(attacks,na.rm=T),sd=sd(attacks,na.rm=T)))
}




#######################
## Coefficient table
#######################

CoefTable <- function(mod,rounding=5){
coefs <- summary(mod)$p.coeff
vcov <- summary(mod)$cov.scaled[1:length(coefs),1:length(coefs)]
alphas <- coefs[1:(length(coefs)/2)]
gammas <- coefs[(length(coefs)/2+1):length(coefs)]
betas <- alphas+gammas
alphas.se <- sqrt(diag(vcov[names(alphas),names(alphas)]))
betas.se <- sqrt(diag(vcov[names(alphas),names(alphas)])+diag(vcov[names(gammas),names(gammas)])+2*diag(vcov[names(alphas),names(gammas)]))


z.alpha <- alphas/alphas.se
z.betas <- betas/betas.se
p.alpha <- 2*pnorm(-abs(z.alpha))
p.betas <- 2*pnorm(-abs(z.betas))
p.alpha <- ifelse(p.alpha<.001,"***",ifelse(p.alpha<.01,"**",ifelse(p.alpha<.05,"*","")))
p.betas <- ifelse(p.betas<.001,"***",ifelse(p.betas<.01,"**",ifelse(p.betas<.05,"*","")))
alphas <- round(alphas,rounding)
alphas.se <- round(alphas.se,rounding)
betas <- round(betas,rounding)
betas.se <- round(betas.se,rounding)
s.edf <- round(summary(mod)$edf,rounding)
s.p <- ifelse(summary(mod)$s.pv<.001,"***",ifelse(summary(mod)$s.pv<.01,"**",ifelse(summary(mod)$s.pv<.05,"*","")))
coef.table <- matrix(NA,nrow=length(alphas)*2+5,ncol=3)
coef.table[1:length(alphas)*2-1,1] <- names(alphas)
coef.table[1:length(alphas)*2,1] <- ""
coef.table[1:length(alphas)*2-1,2] <- alphas
coef.table[1:length(alphas)*2,2] <- paste("(",alphas.se,")",p.alpha,sep="")
coef.table[1:length(alphas)*2-1,3] <- betas
coef.table[1:length(alphas)*2,3] <- paste("(",betas.se,")",p.betas,sep="")
coef.table[(nrow(coef.table)-4):nrow(coef.table)] <- c("Spline","AIC","AUC.in","AUC.out","N")
coef.table[nrow(coef.table)-4,2] <- paste(s.edf,s.p,sep="")
coef.table[nrow(coef.table)-3,2] <- round(mod$aic,2)
coef.table[nrow(coef.table)-2,2] <- round(AUC(mod$y,mod$fitted),4)
coef.table[nrow(coef.table),2] <- summary(mod)$n
colnames(coef.table) <- c("Terms","New","Recurring")
coef.table <- as.data.frame(coef.table)
return(coef.table)
}




#######################
## Coefficient table (FE)
#######################


CoefTableFE <- function(mod,rounding=5){
coefs <- summary(mod)$p.coeff
nFE <- !grepl("OBLAST_NAM",names(coefs))
coefs <- coefs[nFE]
vcov <- summary(mod)$cov.scaled[1:length(coefs),1:length(coefs)]
alphas <- coefs[1:(length(coefs)/2)]
gammas <- coefs[(length(coefs)/2+1):length(coefs)]
betas <- alphas+gammas
alphas.se <- sqrt(diag(vcov[names(alphas),names(alphas)]))
betas.se <- sqrt(diag(vcov[names(alphas),names(alphas)])+diag(vcov[names(gammas),names(gammas)])+2*diag(vcov[names(alphas),names(gammas)]))


z.alpha <- alphas/alphas.se
z.betas <- betas/betas.se
p.alpha <- 2*pnorm(-abs(z.alpha))
p.betas <- 2*pnorm(-abs(z.betas))
p.alpha <- ifelse(p.alpha<.001,"***",ifelse(p.alpha<.01,"**",ifelse(p.alpha<.05,"*","")))
p.betas <- ifelse(p.betas<.001,"***",ifelse(p.betas<.01,"**",ifelse(p.betas<.05,"*","")))
alphas <- round(alphas,rounding)
alphas.se <- round(alphas.se,rounding)
betas <- round(betas,rounding)
betas.se <- round(betas.se,rounding)
s.edf <- round(summary(mod)$edf,rounding)
s.p <- ifelse(summary(mod)$s.pv<.001,"***",ifelse(summary(mod)$s.pv<.01,"**",ifelse(summary(mod)$s.pv<.05,"*","")))
coef.table <- matrix(NA,nrow=length(alphas)*2+5,ncol=3)
coef.table[1:length(alphas)*2-1,1] <- names(alphas)
coef.table[1:length(alphas)*2,1] <- ""
coef.table[1:length(alphas)*2-1,2] <- alphas
coef.table[1:length(alphas)*2,2] <- paste("(",alphas.se,")",p.alpha,sep="")
coef.table[1:length(alphas)*2-1,3] <- betas
coef.table[1:length(alphas)*2,3] <- paste("(",betas.se,")",p.betas,sep="")
coef.table[(nrow(coef.table)-4):nrow(coef.table)] <- c("Spline","AIC","AUC.in","AUC.out","N")
coef.table[nrow(coef.table)-4,2] <- paste(s.edf,s.p,sep="")
coef.table[nrow(coef.table)-3,2] <- round(mod$aic,2)
coef.table[nrow(coef.table)-2,2] <- round(AUC(mod$y,mod$fitted),4)
coef.table[nrow(coef.table),2] <- summary(mod)$n
colnames(coef.table) <- c("Terms","New","Recurring")
coef.table <- as.data.frame(coef.table)
return(coef.table)
}




###############################
## Correlation matrix legend
###############################

legend.v2 <- function (x, y = NULL, legend, fill = NULL, col = par("col"), border = "black", lty, lwd, pch, angle = 45, density = NULL, bty = "o", bg = par("bg"), box.lwd = par("lwd"), box.lty = par("lty"), box.col = par("fg"), pt.bg = NA, cex = 1, pt.cex = cex, pt.lwd = lwd, xjust = 0, yjust = 1, x.intersp = 1, y.intersp = 1, adj = c(0, 0.5), text.width = NULL, text.col = par("col"), merge = do.lines && has.pch, trace = FALSE, plot = TRUE, ncol = 1, horiz = FALSE, title = NULL, inset = 0, xpd, title.col = text.col, title.adj = 0.5,seg.len = 2){
if (missing(legend) && !missing(y) && (is.character(y) || 
is.expression(y))) {
legend <- y
y <- NULL
}
mfill <- !missing(fill) || !missing(density)
if (!missing(xpd)) {
op <- par("xpd")
on.exit(par(xpd = op))
par(xpd = xpd)
}
title <- as.graphicsAnnot(title)
if (length(title) > 1) 
stop("invalid title")
legend <- as.graphicsAnnot(legend)
n.leg <- if (is.call(legend)) 
1
else length(legend)
if (n.leg == 0) 
stop("'legend' is of length 0")
auto <- if (is.character(x)) 
match.arg(x, c("bottomright", "bottom", "bottomleft", 
"left", "topleft", "top", "topright", "right", "center"))
else NA
if (is.na(auto)) {
xy <- xy.coords(x, y)
x <- xy$x
y <- xy$y
nx <- length(x)
if (nx < 1 || nx > 2) 
stop("invalid coordinate lengths")
}
else nx <- 0
xlog <- par("xlog")
ylog <- par("ylog")
rect2 <- function(left, top, dx, dy, density = NULL, angle, 
...) {
r <- left + dx
if (xlog) {
left <- 10^left
r <- 10^r
}
b <- top - dy
if (ylog) {
top <- 10^top
b <- 10^b
}
rect(left, top, r, b, angle = angle, density = density, 
...)
}
segments2 <- function(x1, y1, dx, dy, ...) {
x2 <- x1 + dx
if (xlog) {
x1 <- 10^x1
x2 <- 10^x2
}
y2 <- y1 + dy
if (ylog) {
y1 <- 10^y1
y2 <- 10^y2
}
segments(x1, y1, x2, y2, ...)
}
points2 <- function(x, y, ...) {
if (xlog) 
x <- 10^x
if (ylog) 
y <- 10^y
points(x, y, ...)
}
text2 <- function(x, y, ...) {
if (xlog) 
x <- 10^x
if (ylog) 
y <- 10^y
text(x, y, ...)
}
if (trace) 
catn <- function(...) do.call("cat", c(lapply(list(...), 
formatC), list("\n")))
cin <- par("cin")
Cex <- cex * par("cex")
if (is.null(text.width)) 
text.width <- max(abs(strwidth(legend, units = "user", 
cex = cex)))
else if (!is.numeric(text.width) || text.width < 0) 
stop("'text.width' must be numeric, >= 0")
xc <- Cex * xinch(cin[1L], warn.log = FALSE)
yc <- Cex * yinch(cin[2L], warn.log = FALSE)
if (xc < 0) 
text.width <- -text.width
xchar <- xc
xextra <- 0
yextra <- yc * (y.intersp - 1)
ymax <- yc * max(1, strheight(legend, units = "user", cex = cex)/yc)
ychar <- yextra + ymax
if (trace) 
catn("  xchar=", xchar, "; (yextra,ychar)=", c(yextra, 
ychar))
if (mfill) {
xbox <- xc * 0.8
ybox <- yc * 0.5
dx.fill <- xbox
}
do.lines <- (!missing(lty) && (is.character(lty) || any(lty > 
0))) || !missing(lwd)
n.legpercol <- if (horiz) {
if (ncol != 1) 
warning("horizontal specification overrides: Number of columns := ", 
n.leg)
ncol <- n.leg
1
}
else ceiling(n.leg/ncol)
has.pch <- !missing(pch) && length(pch) > 0
if (do.lines) {
x.off <- if (merge) 
-0.7
else 0
}
else if (merge) 
warning("'merge = TRUE' has no effect when no line segments are drawn")
if (has.pch) {
if (is.character(pch) && !is.na(pch[1L]) && nchar(pch[1L], 
type = "c") > 1) {
if (length(pch) > 1) 
warning("not using pch[2..] since pch[1L] has multiple chars")
np <- nchar(pch[1L], type = "c")
pch <- substr(rep.int(pch[1L], np), 1L:np, 1L:np)
}
}
if (is.na(auto)) {
if (xlog) 
x <- log10(x)
if (ylog) 
y <- log10(y)
}
if (nx == 2) {
x <- sort(x)
y <- sort(y)
left <- x[1L]
top <- y[2L]
w <- diff(x)
h <- diff(y)
w0 <- w/ncol
x <- mean(x)
y <- mean(y)
if (missing(xjust)) 
xjust <- 0.5
if (missing(yjust)) 
yjust <- 0.5
}
else {
h <- (n.legpercol + (!is.null(title))) * ychar + yc
w0 <- text.width + (x.intersp + 1) * xchar
if (mfill) 
w0 <- w0 + dx.fill
if (do.lines) 
w0 <- w0 + (seg.len + +x.off) * xchar
w <- ncol * w0 + 0.5 * xchar
if (!is.null(title) && (abs(tw <- strwidth(title, units = "user", 
cex = cex) + 0.5 * xchar)) > abs(w)) {
xextra <- (tw - w)/2
w <- tw
}
if (is.na(auto)) {
left <- x - xjust * w
top <- y + (1 - yjust) * h
}
else {
usr <- par("usr")
inset <- rep(inset, length.out = 2)
insetx <- inset[1L] * (usr[2L] - usr[1L])
left <- switch(auto, bottomright = , topright = , 
right = usr[2L] - w - insetx, bottomleft = , 
left = , topleft = usr[1L] + insetx, bottom = , 
top = , center = (usr[1L] + usr[2L] - w)/2)
insety <- inset[2L] * (usr[4L] - usr[3L])
top <- switch(auto, bottomright = , bottom = , bottomleft = usr[3L] + 
h + insety, topleft = , top = , topright = usr[4L] - 
insety, left = , right = , center = (usr[3L] + 
usr[4L] + h)/2)
}
}
if (plot && bty != "n") {
if (trace) 
catn("  rect2(", left, ",", top, ", w=", w, ", h=", 
h, ", ...)", sep = "")
rect2(left, top, dx = w, dy = h, col = bg, density = NULL, 
lwd = box.lwd, lty = box.lty, border = box.col)
}
xt <- left + xchar + xextra + (w0 * rep.int(0:(ncol - 1), 
rep.int(n.legpercol, ncol)))[1L:n.leg]
yt <- top - 0.5 * yextra - ymax - (rep.int(1L:n.legpercol, 
ncol)[1L:n.leg] - 1 + (!is.null(title))) * ychar
if (mfill) {
if (plot) {
fill <- rep(fill, length.out = n.leg)
rect2(left = xt+.22, top = yt + ybox/2, dx = xbox, dy = ybox, 
col = fill, density = density, angle = angle, 
border = border)
}
xt <- xt + dx.fill
}
if (plot && (has.pch || do.lines)) 
col <- rep(col, length.out = n.leg)
if (missing(lwd)) 
lwd <- par("lwd")
if (do.lines) {
if (missing(lty)) 
lty <- 1
lty <- rep(lty, length.out = n.leg)
lwd <- rep(lwd, length.out = n.leg)
ok.l <- !is.na(lty) & (is.character(lty) | lty > 0)
if (trace) 
catn("  segments2(", xt[ok.l] + x.off * xchar, ",", 
yt[ok.l], ", dx=", seg.len * xchar, ", dy=0, ...)")
if (plot) 
segments2(xt[ok.l] + x.off * xchar, yt[ok.l], dx = seg.len * 
xchar, dy = 0, lty = lty[ok.l], lwd = lwd[ok.l], 
col = col[ok.l])
xt <- xt + (seg.len + x.off) * xchar
}
if (has.pch) {
pch <- rep(pch, length.out = n.leg)
pt.bg <- rep(pt.bg, length.out = n.leg)
pt.cex <- rep(pt.cex, length.out = n.leg)
pt.lwd <- rep(pt.lwd, length.out = n.leg)
ok <- !is.na(pch) & (is.character(pch) | pch >= 0)
x1 <- (if (merge && do.lines) 
xt - (seg.len/2) * xchar
else xt)[ok]
y1 <- yt[ok]
if (trace) 
catn("  points2(", x1, ",", y1, ", pch=", pch[ok], 
", ...)")
if (plot) 
points2(x1, y1, pch = pch[ok], col = col[ok], cex = pt.cex[ok], 
bg = pt.bg[ok], lwd = pt.lwd[ok])
}
xt <- xt + x.intersp * xchar
if (plot) {
if (!is.null(title)) 
text2(left + w * title.adj, top - ymax, labels = title, 
adj = c(title.adj, 0), cex = cex, col = title.col)
text2(xt, yt, labels = legend, adj = adj, cex = cex, 
col = text.col)
}
invisible(list(rect = list(w = w, h = h, left = left, top = top), 
text = list(x = xt, y = yt)))
}


#######################################
## Outreg function (by Paul Johnson)
#######################################


outreg <- function(incoming, title="My Regression", label="", modelLabels=NULL, varLabels=NULL, tight=TRUE, showAIC=TRUE, lyx=TRUE){
modelList <- NULL
## was input just one model, or a list of models? ###
if ( "lm" %in% class(incoming)) { ##just one model input
nmodels <- 1
modelList <- list(modl1=incoming)
} else {
nmodels <- length(incoming)
modelList <- incoming
}

##TODO modelLabels MUST have same number of items as "incoming"

## Get a regression summary object for each fitted model
summaryList <- list()
fixnames <- vector()
myModelClass <- vector()
i <- 1
for (model in modelList){
summaryList[[i]] <- summary(model)
fixnames <- unique( c( fixnames, names(coef(model))))
myModelClass[i] <- class(model)[1]
i <- i+1
}


###If you are just using LaTeX, you need these

if (lyx == FALSE){
cat("\\begin{table}\n ")
cat("\\caption{",title,"}\\label{",label,"}\n ")
}
cat("\\begin{center}\n ")
nColumns <- ifelse(tight, 1+nmodels, 1 + 2*nmodels)
cat(paste("\\begin{tabular}{*{",nColumns,"}{l}}\n ", sep=""))
cat("\\hline\n ")

### Put model labels on top of each model column, if modelLabels were given

if (!is.null(modelLabels)){
cat(" ")
for (modelLabel in modelLabels){
if (tight == T) {
cat(paste("&", modelLabel))
}else{
cat(paste("&\\multicolumn{2}{c}{",modelLabel,"}",sep=""))
}
}
cat (" \\\\\n ")
}

### Print the headers "Estimate" and "(S.E.)", output depends on tight or other format

if (tight == T){
cat(" ")
for (i in 1:nmodels) { cat (" & Estimate ") }
cat(" \\\\\n")
cat(" ")
for (i in 1:nmodels) { cat (" & (S.E.) ") }
cat(" \\\\\n")
}else{
cat(" ")
for (i in 1:nmodels) { cat (" & Estimate & S.E.") }
cat(" \\\\\n")
}
cat("\\hline \n \\hline\n ")

### Here come the regression coefficients

for (regname in fixnames){
if ( !is.null(varLabels[[regname]]) ) { cat(paste("",varLabels[[regname]]), sep="")}
else {cat(paste("", regname), sep="")}
for (model in modelList) {
est <- coef(model)[regname]
se <- sqrt(diag(vcov(model)))[regname]
if ( !is.na(est) ) {
cat (paste(" & ", round(est,5)))
pval <- pt(abs(est/se), lower.tail=F, df = model$df.residual)
if(is.na(pval)==0){
if (pval < 0.001) cat("***")
if (pval >= 0.001 & pval < 0.01) cat("**")
if (pval >= 0.01 & pval < 0.05) cat("*")}
if (tight == F) {
cat (paste(" & (", round(se,5),")",sep=""))
}
} else {
cat (" & . ")
if (tight == F) cat (" & " )
}

}
cat (" \\\\\n ")
if (tight == T){
for (model in modelList) {
est <- coef(model)[regname]
if (!is.na(est)) cat (paste(" & (",round(sqrt(diag(vcov(model)))[regname],4)),")",sep="")
else cat(" & ")
}
cat (" \\\\\n ")
}
}
cat("\\hline \n")

### Print a row for the number of cases

cat(paste("N"), sep="")
for (model in summaryList) {
myDF <- sum( model$df[-3] ) #omit third value from df vector
cat (paste(" & ", myDF))
if (tight == F) cat(" &")
}
cat (" \\\\\n ")

### Print a row for the root mean square error

if ("lm" %in% myModelClass) {
cat(paste("$RMSE$"),sep="")
for (model in summaryList) {
cat( paste(" &", if(is.numeric(model$sigma)) round(model$sigma,4)))
if (tight == F) cat(" &")
}
cat (" \\\\\n ")
}

### Print a row for the R-square

if ("lm" %in% myModelClass) {
cat(paste("$R2$"),sep="")
for (model in summaryList) {
cat( paste(" &", if(is.numeric(model$r.square))round(model$r.square,4)))
if (tight == F) cat(" &")
}
cat (" \\\\\n ")
}

## Print a row for the model residual deviance

if ("glm" %in% myModelClass) {
cat(paste("$Deviance$"),sep="")
for (model in summaryList) {
cat (paste(" &", if(is.numeric(model$deviance))round(model$deviance,4)))
if (tight == F) cat(" &")
}
cat (" \\\\\n ")
}

### Print a row for the model's fit, as -2LLR

if ("glm" %in% myModelClass) {
cat (paste("$-2LLR (Model \\chi2)$"),sep="")
for (model in modelList) {
if (is.numeric(model$deviance)){
n2llr <- model$null.deviance - model$deviance
cat (paste(" &", round(n2llr,4)))
gmdf <- model$df.null - model$df.residual + 1
if (pchisq(n2llr, df= gmdf, lower.tail=F) < 0.001) {cat ("***")}
if (pchisq(n2llr, df= gmdf, lower.tail=F) >= 0.001 & pchisq(n2llr, df= gmdf, lower.tail=F) < 0.01) {cat ("*")}
if (pchisq(n2llr, df= gmdf, lower.tail=F) >= 0.01 & pchisq(n2llr, df= gmdf, lower.tail=F) < 0.05) {cat ("*")}
}
else {
cat (" &")
}
if (tight == F) cat(" &")
}
cat (" \\\\\n ")
}


## Print a row for the model's fit, as -2 LLR
### Can't remember why I was multiplying by -2

if (showAIC == T) {
cat(paste("$AIC$"),sep="")
for (model in modelList) {
cat (paste(" &", if(is.numeric(AIC(model)))round(AIC(model),3)))
if (tight == F) cat(" &")
}
cat (" \\\\\n ")
}
cat("\\hline\\hline\n")
cat ("* $p \\le 0.05$, ** $p \\le 0.01$, *** $p \\le 0.001$")
cat("\\end{tabular}\n")
cat("\\end{center}\n")
if (lyx == FALSE){
cat("\\end{table}\n")
}
}
